/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of UHM (Unassembled HyperMatrix) sparse direct solver.
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/
#include "uhm.hxx"

#define UHM_ERROR_TOL 1.0e-5

int main (int argc, char **argv)
{

  // input check
  if (argc != 6) {
    printf("Try : %s [n_thread][n_parts][decomposition][blocksize][input_file]\n\n",
           argv[0]);
    printf(" ==== decomposition method ====\n\n");
    printf(" ** Symmetric Factorization **\n");
    printf("    UHM_CHOL      ( 1),   UHM_CHOL_PIV ( 2),   UHM_CHOL_INCPIV ( 3) \n");
    printf("    UHM_LDL_NOPIV ( 4),   UHM_LDL_PIV  ( 5),   UHM_LDL_INCPIV  ( 6) \n");
    printf("\n");
    printf(" ** General Structure Symmetric Factorization **\n");
    printf("    UHM_LU_NOPIV  ( 7),   UHM_LU_PIV   ( 8),   UHM_LU_INCPIV   ( 9) \n");
    printf("    UHM_QR        (10) \n");
    return 0;
  }

  int n_threads, n_parts, decomposition, blocksize;
  char *filename;

  n_threads     = atoi( (argv[1]) );
  n_parts       = atoi( (argv[2]) );
  decomposition = atoi( (argv[3]) );
  blocksize     = atoi( (argv[4]) );
  filename      = argv[5];

  double t_base, t_part_mesh, t_build_tree, t_decompose;

  FLA_Init();

  uhm::set_num_threads    (n_threads);
  uhm::set_hier_blocksize (blocksize);

  int sym      = ( decomposition < UHM_LU_NOPIV ? UHM_SYMMETRY : UHM_UNSYMMETRY);
  int datatype = UHM_REAL, n_rhs = 1, is_schur = false;

  uhm::Net  net = new uhm::Net_;
  uhm::Mesh m   = net->get_self();

  printf( "BEGIN : Import mesh from file : %s\n", filename );
  m->import_file( filename );
  printf( "END   : Import \n" );

  printf( "BEGIN : Partition mesh \n" );
  t_base       = uhm::timer();
  uhm::partition_mesh(m, n_parts);
  t_part_mesh  = uhm::timer() - t_base;
  printf( "END   : Partition mesh %lf\n", t_part_mesh );

  m->lock();
  net->lock();

  for (int i=0;i<net->get_n_children();++i) {
    uhm::Mesh c = net->get_child(i);

    printf( "BEGIN : Build tree of child %d mesh \n", i );
    t_base       = uhm::timer();
    uhm::build_tree(c);
    t_build_tree = uhm::timer() - t_base;
    printf( "END   : Build tree of child %d mesh %lf\n", i, t_build_tree );
    
    c->set_symmetry(sym);
    c->lock();

    printf(" n_nodes %d, n_elts %d\n", c->get_n_nodes(), c->get_n_elements());

    c->create_matrix_without_buffer( datatype, n_rhs );
    c->create_matrix_buffer(is_schur);
    if (decomposition < UHM_LU_NOPIV) 
      c->random_spd_matrix();
    else 
      c->random_matrix();
    c->set_rhs();
  }

  printf("BEGIN : Decomposition of children mesh with decomposition %d\n",
         decomposition);

  t_base      = uhm::timer();
  net->decompose_local_with_free(decomposition);
  t_decompose = uhm::timer() - t_base;
  printf("END   : Decomposition of children mesh with decomposition %d, %lf\n", 
         decomposition, t_decompose);


  //net->disp();

  delete net;

  FLA_Finalize();

  return 0;
}


